function z = SS_eq( alpha1,alpha2,alpha3,z_grid,lambda_grid,im,ih )
%UNTITLED2 Summary of this function goes here
%   Detailed explanation goes here
if (im>(alpha1+alpha3)*lambda_grid(1))&&(ih>(alpha2+alpha3)*lambda_grid(1))
    z=[0,0];
else
    zm_can = interp1((alpha1+alpha3)*lambda_grid,z_grid,im);
    zh_can = interp1((alpha2+alpha3)*lambda_grid,z_grid,ih);
    Ah = alpha2*lambda_grid(1)+alpha3*interp1(z_grid,lambda_grid,zm_can);
    Am = alpha1*lambda_grid(1)+alpha3*interp1(z_grid,lambda_grid,zh_can);
    
    if Ah<ih
        z = [zm_can,0];
    elseif Am<im
        z = [0,zh_can];
    else
        AA1 = im-ih+alpha2*lambda_grid;
        zm_zh = interp1(alpha1*lambda_grid,z_grid,AA1);
        ID = ~isnan(zm_zh);
        zm_zh = zm_zh(ID);
        z_grid1 = z_grid(ID);
        AA2 = alpha2*lambda_grid(ID)+alpha3*interp1(z_grid,lambda_grid,min(z_grid1+zm_zh,max(z_grid)))-ih;
        [~,ID] = min(abs(AA2));
        %zh = fzero(@(x)interp1(z_grid1,AA2,x),[min(z_grid1),max(z_grid1)]);
        zh =z_grid1(ID); %fmincon(@(x)interp1(z_grid1,AA2,x)^2,zm_zh(ID),[],[],[],[],min(z_grid1),max(z_grid1));
        z = [interp1(z_grid1,zm_zh,zh),zh];
        
    end
    
end


end

